subroutine symspec(kb,k,bgx,bgy,ssz5,nkhz)

!  kb    tells which buffer (0 or 1) is in use
!  k     pointer to the most recent new data
!  px    power in x channel
!  py    power in y channel
!  s     spectrum
!  nkhz  kHz of center frequency, e.g. 125 if fcenter=144.125

  parameter (NSMAX=60*96000)          !Total sample intervals per minute
  parameter (NFFT=32768)              !Length of FFTs
  real ssz5(NFFT)
  real*8 hsym
  complex cx(NFFT),cy(NFFT)
  complex z
  include 'datcom.f90'
  common/gcom0/nqsofreq
  data ncall/0/
  save ncall

  hsym=2048.d0*96000.d0/11025.d0      !Samples per JT65 half-symbol
  nhsym=hsym
  jb=k                                !End index 
  if(jb.gt.NSMAX) jb=jb-NSMAX
  ja=jb-nhsym+1                       !Start index
  if(ja.lt.1) ja=1

  px=0.
  py=0.
  i=0
  do j=ja,jb                          !Copy data into cx, cy
     x1=dd(1,j,kb)
     x2=dd(2,j,kb)
     x3=dd(3,j,kb)
     x4=dd(4,j,kb)
     i=i+1
     cx(i)=cmplx(x1,x2)
     cy(i)=cmplx(x3,x4)
     px=px + x1*x1 * x2*x2            !Compute power
     py=py + x3*x3 * x4*x4
  enddo
  nsum=jb-ja+1
  rmsx=sqrt(0.5*px/nsum)
  rmsy=sqrt(0.5*py/nsum)
  bgx=0.
  bgy=0.
  if(rmsx.gt.1.0) bgx=log10(rmsx)/log10(32768.0)
  if(rmsy.gt.1.0) bgy=log10(rmsy)/log10(32768.0)
  if(bgx.gt.1.0) bgx=1.0
  if(bgy.gt.1.0) bgy=1.0

  u=0.1
  nclimit=nint(2.0/u)
  pwr=(rmsx/100.0)**2 + (rmsy/100.0)**2
  if(ncall.lt.1000000) ncall=ncall+1
  if(ncall.gt.nclimit .and. pwr.gt.3.0*avgpwr) then
!     write(84,3001) ncall,k/96000.0,pwr,avgpwr
!3001 format(i10,f10.3,2f10.1)
!     call flush(84)
     do j=ja,jb                          !Blank the bad data
        dd(1:4,j,kb)=0.
     enddo  
     cx(:i)=0.
     cy(:i)=0.
  endif
  up=pwr-avgpwr
  if(ncall.gt.nclimit) up=min(up,0.2*avgpwr)
  avgpwr=avgpwr + u*up

  cx(i+1:)=0.                         !Pad with zeros
  cy(i+1:)=0.
  call four2a(cx,NFFT,1,1,1)          !Compute FFTs
  call four2a(cy,NFFT,1,1,1)

  do i=1,NFFT
     sx=real(cx(i))**2 + aimag(cx(i))**2        
     sy=real(cy(i))**2 + aimag(cy(i))**2
     z=cx(i)*conjg(cy(i))
     q=sx - sy
     u=2.0*real(z)
     ssz5(i)=0.707*sqrt(q*q + u*u)    !Spectrum of linear polarization
  enddo

  nkhz=1000*(fcenter-int(fcenter))

  return
end subroutine symspec
